Sheared Ising models in three dimensions 



Alfred Hucht and Sebastian Angst 



O 



43 

o 



Fakultat far Physik, Universitdt Duisburg-Essen, D-47048 Duisburg 



PACS 05.70.Ln- 
PACS 68 . 35 . Af - 
PACS 05.50.+q- 

Abstract -The nonequilibrium phase transition in sheared three-dimensional Ising models is in- 
vestigated using Monte Carlo simulations in two different geometries corresponding to different 
shear normals. We demonstrate that in the high shear limit both systems undergo a strongly 
anisotropic phase transition at exactly known critical temperatures T c which depend on the di- 
rection of the shear normal. Using dimensional analysis, we determine the anisotropy exponent 
= 2 as well as the correlation length exponents i/y = 1 and z/j_ = 1/2. These results are verified 
by simulations, though considerable corrections to scaling are found. The correlation functions 
perpendicular to the shear direction can be calculated exactly and show Ornstein-Zernike behavior. 



a 
o 
o 



> 
o 
r-- 

On 

m 

• 

o 



Introduction. — While the occurrence of nonequilib- 
rium phase transitions is ubiquitous in nature, its investi- 
gation in the framework of nonequilibrium statistical me- 
chanics is intricate and restricted to a few simple mod- 
els, like the driven lattice gas (DLG) y}j3] or, recently, 
to the driven two-dimensional Ising model |4j. In this 
model the system is cut into two halves parallel to one 
axis and moved along this cut with the velocity v. The 
model exhibits energy dissipation and subsequently fric- 
tion due to spin correlations, which also occurs in a suit- 
able Heisenberg model [5}j8] and, of interest for the current 
context, undergoes a nonequilibrium (surface) phase tran- 
sition. The latter has been investigated analytically and 
with Monte Carlo (MC) simulations for various geome- 
tries 1 9 1. Since then, this model has been generalized to 
the driven Potts models JlO) , and finite-size effects were 
calculated analytically in the driven Ising chain [TT] , 

A lot of similarities and comparable critical behavior 
between the Ising model with friction and the very fa- 
mous and well investigated DLG have been found Jl2] , 
Both models are characterized by a critical temperature 
T c , which increases with the driving strength, the field 
and the shift or shear velocity v, respectively, and satu- 
rates in the high driving limit. For diverse geometries of 
the Ising model with friction, the critical temperature has 
been calculated analytically for v — >• oo |9|. 

Moreover it was discovered that the DLG and two- 
dimensional sheared Ising systems with non-conserved or- 
der parameter |12]jl4| show strongly anisotropic critical 
behavior, with direction dependent correlation length ex- 



ponents v\\ and For the 2d and 1+ld geometry of 
the Ising model with shear the same exponents z/y =3/2 
and v±_ = 1/2 |12 as in the two-dimensional DLG have 
been determined. Additionally finite velocities v have been 
studied and it was found that for all finite v the 2d and 
1+ld model cross-over from isotropic Ising like behavior 
to strongly anisotropic mean-field behavior in the thermo- 
dynamic limit, demonstrating that the external drive is a 
relevant perturbation. 

In the following we extend the investigations to three- 
dimensional models in two different shear geometries and 
focus on the high shear velocity limit v — »• oo. Both repre- 
sent three-dimensional sheared models and they are there- 
fore experimentally accessible in the framework of sheared 
binary liquids [l5}{l8] , albeit the order parameter is not 
conserved here. Using dimensional analysis, we predict 
the correlation length exponents for arbitrary dimension 
d. These predictions are verified by simulations, however 
we find strong corrections to scaling at small system sizes. 

Model. — The systems considered in this work are 
denoted 2+ld and l+2d and are shown in Fig. [I] for a 
classification see Ref. |9|. In the 2+ld geometry shear is 
applied such that two-dimensional Ising models are moved 
relative to their upper (lower) neighboring layer with ve- 
locity v (—v) along an axis. In the following we denote 
the direction parallel to the shear with ||, the direction 
perpendicular to the planes with J_i and the inplane di- 
rection perpendicular to the shear direction with _L 2 - The 
model contains L± 1 x Lj_ 2 x Ln spins (lattice sites), where 
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Figure 1: Sketches of the systems considered in this work. On 
the left hand side the 2+ld system and on the right hand side 
the l+2d system is shown. The gray regions represent the 
magnetic systems and the green (dark) regions are the mov- 
ing boundaries, while the arrows indicate the motion of the 
subsystems. 



The corresponding Hamiltonian reads 



pH(t) 



o~kim(^wak 7 i,m+i + 

k=l 1=1 rn=l 
~ ^_L[o"fe,Z+l,m+A(t) + 0"fc+l,Z,m+A(t)]) j ( 4 ) 

=: K±. For v — » oo we set / = 4 and 
e 2K c,\\ from the one-dimensional Ising 
to get, for J\\ = J± = 1, the critical 



where K± 1 = K±, 

USe Xeq(^c,||) = 

model in Eq. Q 
temperature 



we choose L± t = L± 2 =: L± throughout this work, and 
periodic boundary conditions are applied in all directions. 
The shear velocity v corresponds to a shear rate, which is 
often denoted as 7 [l3{[l4] . Using the notation (±1 J_2 ||) 
for directions, the shear is in (OOl)-direction and the shear 
normal is in (lOO)-direction. 

A finite shear velocity v is implemented by shifting 
neighboring layers v times by one lattice constant dur- 
ing one MC step (for details see |4,9|). A simplification of 
the implementation is yielded by reordering the couplings 
between moved layers instead, and by introducing a time- 
dependent displacement A(t) = vt we get the Hamiltonian 



l+2d 



1 - 1 - z II 

/3H(t) = ^2 0klrn[K\\<7k,l,m+l + 

k=l 1=1 rn=l 

where = j3 J M is the reduced nearest neighbor coupling 
with \i = {J_i,J_2, ||}, and f3 = In the following 

we concentrate on the infinite shear velocity limit v —> 00, 
which can easily be implemented by choosing 1 < A(t) < 
L|| randomly. In this limit an analytical calculation |9] 
yield the equation 



\)f tanhif. 



(2) 



from which we can determine the critical temperature, 
where Xeq is the zero field equilibrium susceptibility of 
the subsystems moved relative to each other and / the 
number of fluctuating adjacent fields. Here Xeq of the 
two-dimensional Ising model is required, which has been 
calculated to higher than 2000 th order by an polynomial 
algorithm |19 . Using / = 2 and Jm = Jj_ x = J± 2 = 1 we 
get 



T c 2+ld (oo) = 5.2647504145147435505980. 



(3) 



The second considered geometry l+2d is similar to the 
previous case, but now the shear normal is in the (In- 
direction. As a consequence, all four perpendicular cou- 
pling partners of a spin a are in neighboring shear planes. 



(00) 



ln[§(5 + V4l)] 



5.642611138.... 



(5) 



which notably is different from Eq. Hence the critical 
temperature depends on the direction of the shear normal. 

In MC simulations of nonequilibrium models the critical 
temperature often depends on the used acceptance rates 
1 20 1 . It has been shown that the multiplicative rate |9 



Pm P (AE) = e 



}(AE-E mill ) 



(6) 



with the energy change AE and the minimal energy 
change AE m { n = min { AE} must be used in order to re- 
produce the critical temperatures Eqs. (|3|[5|. 



Anisotropic scaling. — Our aim is to proof that both 
models exhibit a strongly anisotropic phase transition and 
calculate the corresponding exponents. Such a phase tran- 
sition is characterized by bulk correlation lengths £ M di- 
verging with direction dependent critical exponents v ^ at 
criticality Q 

^(t)%°i^, (?) 

with direction /i = {_Li, _L 2 > ||}> amplitude £ M , and reduced 
critical temperature t — T/T c — 1. Usually one defines 
the anisotropy exponent 6 = v\\/v±, which is 6 = 1 for 
isotropic scaling and 7^ 1 for strongly anisotropic scal- 
ing [2{ |2l}|24] . As mentioned above, the phase transitions 
of the Ising model with friction in the 2d and the 1+ld 
geometry become strongly anisotropic for v > in the 
thermodynamic limit, with = 3 [12] . 

In Ref. 1 12 it was shown that the application of a stripe 
geometry L± — >• 00 with finite Ln is an appropriate way 
to determine the anisotropy exponent and subsequently 
the correlation length exponents. Hence we measure the 
perpendicular correlation function 



G±(L\\]r±) = (croooov 



(8) 



at the critical point T c , from which we can determine the 
correlation lengths £ M with \i = {J_i,J_2} as shown below 
(in the following the index \i only represents the perpen- 
dicular directions ±1 and ±2)- Note that by symmetry 



1 Throughout this work the symbol ^ 
cally equal" in the respective limit, e.g., 
limz^oo f(L)/g(L) = 1 . 



means 
f(L) 



"asymptoti- 

- g{L) «*■ 
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G±(Ln,r± 1 ) = G±(Lu , r_L 2 ) for the l+2d system. From 
£ M we can then determine using the relation |23}|25| 



1/9 



(9) 



The above-mentioned stripe geometry is a film geometry 
in three dimensions, and we choose L±/£±(L\\) > 10 suf- 
ficient for our purpose |12 1. 

Dimensional analysis. — For v — >• oo it was shown 
in Ref. [9]| that the 1+ld model can be mapped onto an 
equilibrium system consisting of one-dimensional chains 
that only couple via fluctuating magnetic fields. Due to 
the stripe geometry with short length Ln and the peri- 
odic boundary conditions in parallel direction the magne- 
tization m(x) with x = (x_l,#||) is homogeneous in this 
direction, and parallel correlations are irrelevant. Hence 
we can use the zero mode approximation in this direction, 
leading to an order parameter m = m(xj_) only. 

The resulting Ginzburg-Landau- Wilson (GLW) Hamil- 
tonian 



m = h 



dx 



4-1 



to 1 /t— 7 \ o ^ i 

-m + - Vm + — t m/ 
2 2 4! 



(10) 



can, however, not be mapped onto a Schrodinger equation 
for systems with d > 2 as done in Ref. |12 , as the (d—1)- 
dimensional integral cannot be interpreted as a time in- 
tegral. Instead we use dimensional analysis in order to 
predict the critical exponents: starting from the GLW 
Hamiltonian (10) in d dimensions we eliminate Ln with 



the substitution 



rh L\ 



-l/(5-d) 



X_L 



xL 



l/(5-d) 



t -> tL\ 



-2/(5-d) 



to get the (d— l)-dimensional Hamiltonian 



J dx^ 1 



t n 1 

-fa 1 + - 



(Vm) 



4! 



~ 4 

m 



(lla) 
(lib) 
(11c) 

(12) 



with rh = m(x). From Eqs. ([TT)3,c) we directly read off 
the exponents 



5 — d, 



5-d 



1 

2' 



(13) 



reproducing the results for <i = 1 |9 and d = 2 |12 and 



, 12l 



fulfilling the generalized hyperscaling relation |26 

+ (d - 1)i/_l = 2 - a (14) 
with a = (9II12I. For our case d = 3 we find 



1, 



(15) 



while for d > 4 we predict isotropic or weakly anisotropic 
behavior with = 1 and v\\ = v± = 1/2, as then the upper 
critical dimension d c = 4 is reached and the shear becomes 
an irrelevant perturbation. 



Correlation functions. 



The perpendicular corre- 



lation function can be calculated from Eq. (12) using a 
Gaussian approximation, which is valid, since we investi- 
gate the system at the critical temperature of the bulk, 
which is higher than the the critical temperature of the 
studied film geometry. Setting u = in Eq. (12) and us- 



ing £ oc t 1//2 we get the Ornstein-Zernike structure factor 

1 



S(k) 



~k 2 + 1 



(16) 



In our case the dimension is d = 3, and a Fourier trans- 
formation yields the correlation function 



Using G ex rh 2 oc L 



G(r)<xK (r/£). 
-i/i/ 



(17) 

and back-substituting with 



Eqs. (11) gives the result 
G(Lii;rx) oc L~ 



K [r ± /a(h 



(18) 



for the perpendicular correlation function of the GLW 
Hamiltonian (10), with modified Bessel function of the 



second kind Kq. 

The 2+ Id geometry is weakly anisotropic in perpendic- 
ular direction at least for different couplings J± 1 ^ J^ 2 , 
i.e., the correlation lengths £ > ± 1 and £u have same expo- 
nent v±_ but different amplitudes £ M [23j. This anisotropy 
can be removed by the rescaling 



(19) 



with amplitude from Eq. (PI). Now the perpendicular 



directions are isotropic and we can use Eq. (18) to get the 
final result 



G ± {L r ,rr)*GL» 1/v *K [rML n )] 



(20) 



for the two directions ii = J_i and ±2- Here we already 
have back-substituted with Eq. (19). Note that especially 
in the 2+ Id case the amplitude G should not depend on 
the direction \i. 

Results. — We measured G±_{L\\\r^) at criticality for 
both models using extensive Monte Carlo simulations and 
fitted the results against Eq. po] ) to get shown in 

Fig. [2] As in the 1+ld case we find corrections to scal- 
ing for L || ;$ 300 which are problematic in these three- 
dimensional cases as we cannot simulate systems larger 



Model 






G 


co 


l+2d 
2+ld 


_L 
J-i 


0.254(5) 
0.320(5) 
0.331(5) 


0.93(1) 
0.85(1) 
0.85(1) 


14.(1) 
12.(1) 
12.(1) 



Table 1: Amplitudes and corrections to scaling parameter Co 
for both models. 
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Figure 2: Perpendicular correlation lengths ^(Ly) for the 
l+2d geometry (red circles), the 2+ld geometry in the -Li- 
direction (green diamonds) and in the ±2-direction (blue 
squares) at criticality. The statistical error is smaller than the 
symbol size. Due to corrections to scaling small systems have 
effective anisotropy exponent # e ff ~ 3 (dotted line), which is 
obtained from the logarithmic derivative and shown exemplary 
for system l+2d in the inset. 



than L|i = 1024. Hence we have to introduce a lattice cor- 
rection term in the perpendicular correlation length and 
improve relation ([9| using the ansatz 



Figure 3: Rescaled correlation function G±(L\\ 



with fi 



Uh)=ML\\+coLl /2 + ...) 



si/e 



(21) 



with = 2, which gives the best fit to the data. From 
the numerical data we find the amplitudes and G as 
well as the correction parameter cq listed in Tab. [T] and 
the resulting fit is shown as solid line in Fig. [2] For large 
systems the curve approaches the theoretical limit Eq. ([9| 
with slope 0~ x = 1/2. Note that for small Ly < 64 we 
could also find a reasonable data collapse with exponent 
#eff = 3 (dotted line). 

The resulting rescaled correlation functions for both 
models are presented in Fig. [3] In all cases the y-axis can 
be rescaled with L\\ as predicted, without notable correc- 
tions. We find a convincing data collapse onto the mean- 
field correlation function Kq(t/£) from Eq. (20). For small 
distances r_\_ 2 = 0(1) the correlation function G±(L\\] r± 2 ) 
differs from Eq. p0| ) due to the inplane nearest neighbor 
interactions. 

Now we comment on the four-dimensional geometry 
l+3d, with decouples to a three-dimensional array of in- 
teracting chains, with / = 6 in Eq. We performed test 
simulations for system sizes up to 32^ x 32 and found very 
strong, possibly logarithmic corrections to scaling. From 
the scaling behavior of the available data we estimate that 
system sizes L\\,L± > 1000 would be required to find the 
correct scaling behavior. 

Finally, we extend the dimensional analysis to the gen- 
eral case of a d-dimensional hyper-cubic sheared lattice 
with d\\ driven dimensions and d± perpendicular dimen- 
sions. We again must distinguish between the d± 1 dimen- 



{J_,J_i,J_2} for both models at criticality. We show varying 
system extensions L\\ {8,16,32,64,128,256,512,1024} for both 
cases. A rescaling of the x-axis with ^(Ly) and of the y-axis 
with L\\ results in an excellent data collapse, verifying = 2 
and z/ii = 1. The solid lines represent the calculated Ornstein- 



Zernike correlation function, Eq. (20). Note that we multiplied 



the collapsed data with different factors as indicated in order 
to show them in one plot. 



sions normal to the shear and d± 2 "inplane" dimensions 
without shear motion, with d± = d± 1 + d± 2 . The criti- 
cal temperature T c at infinite shear velocity v is given by 
Eq. (J2|, with the equilibrium zero field susceptibility Xeq 
of the deq-dimensional system having / fluctuating fields 
at each lattice point, where d eq = d\\ + d± 2 , and / = 2d± 1 . 



From a simple generalization of Eq. (13) we find the ex- 
ponents 



4-dj 
d\\ 



4-d A 
2dn 



1 

r 



(22) 



fulfilling the hyperscaling relation d\\V\\ + d±v± = 2. 

We conclude with a tabular summary of the found ex- 
ponents and critical temperatures T c at infinite driving 
velocity v given in Table [2| including two cases denoted 
"mix" where we assumed a suitable two-dimensional mo- 
tion of the interacting planes. These systems have d\\ =2, 
but notwithstanding the same T c as the corresponding sys- 
tems with unidirectional motion at infinite v. For the lay- 
ered case 2+ld m we predict the exponents = 3/2 and 
v\\ = 3/4. A test of these predictions is left for future 
work. 

Conclusion. — We investigated the phase transition 
of three-dimensional Ising models with shear and two dif- 
ferent shear normals by means of Monte Carlo simula- 
tions. In the limit of infinitely high shear velocity v we 
found a critical temperature T c (oo) that depends on the 
direction of the shear normal. At criticality, strongly 
anisotropic diverging correlation lengths, with exponents 
uu = 1 and v± = 1/2 occur, leading to an anisotropy 
exponent = 2, which confirms the results of a dimen- 
sional analysis of the corresponding Ginzburg-Landau- 
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model 


d 








d ±2 


e i/,, 


/ 




T c (oo)/J 




id 


1 


1 


- 


- 


- 


- 


2 


1 


1 


2.2691853. . . 


-a 


2d 


2 


1 


1 





1 


3 


3/2 


1 


2 


4.0587824. . . 


> 


3d 


3 


1 


2 





2 


2 


1 


1 


3 


5.983835(1) 


a 


2d b 


1 


1 


- 


- 


- 


- 


2 


1 


2 


2.6614725. . . 




3d b 


2 


1 


1 





1 


3 


3/2 


1 


3 


4.8(1) 




1+ld 


2 


1 


1 


1 





3 


3/2 


2 


1 


3.4659074. . . 


heare 


2+ld 


3 




2 


1 


1 


2 


1 


2 


2 


5.2647504. . . 


l+2d 


3 




2 


2 





2 


1 


4 


1 


5.6426111... 


CO 


l+3d 


4 




3 


3 





1 


1/2 


6 


1 


7.728921. . . 


X 


2d m 


2 


2 










1 


1 


1 


4.0587824. . . 


a 


2+ld m 


3 


2 


1 


1 





3/2 


3/4 


2 


2 


5.2647504. . . 



Table 2: Relevant dimensions, exponents and parameters of the considered models, as defined in the text. For a classification 

see |9j. 



Wilson Hamiltonian. Furthermore, the dimensional anal- 
ysis captures the anisotropy exponents as well as the cor- 
relation length exponents of the previously studied two- 
dimensional cases 1 12 and the parallel correlation length 
exponent of the one-dimensional cases [9]. Predictions for 
two-dimensional shear directions also result from the di- 
mensional analysis, leading to the exponents = 3/2 and 
z/M = 3/4 in a three-dimensional model. Fluctuations per- 
pendicular to the shear were shown to be Gaussian, result- 
ing in a correlation function with Ornstein-Zernike behav- 
ior. Additionally, in the case of the 2+ Id geometry we 
found weakly anisotropic perpendicular correlations. As 
for v = the 2+ Id and the l+2d geometry reduce to 
the three-dimensional equilibrium Ising model, we expect 
a cross-over from this case to strongly anisotropic mean- 
field behavior similar to the 1+ld geometry. In Ref. [12] 
an expensive analysis for finite velocities has been done 
leading to a crossover scaling, pointing out that all v ^ 
provoke strongly anisotropic mean- field behavior, which is 
expected to occur in the current systems as well. How- 
ever, we did not proof this in detail, due to the additional 
complexity in three-dimensional systems. 
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